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We study the non-reversibility of molecular dynamics trajectories arising from the amplification of rounding 
errors. We analyse the causes of such behaviour and give arguments, indicating that this does not pose a significant 
problem for Hybrid Monte Carlo computations. We present data for pure 517(3) gauge theory and for QCD with 
dynamical fermions on small lattices to illustrate and to support some of our ideas. 



The theory of the Hybrid Monte Carlo (HMC) 
algorithm assumes the exact reversibility of its 
molecular dynamics (MD) trajectories. Leapfrog 
integration guarantees this unless the initial con- 
jugate gradient (CG) vector is chosen in time 
asymmetric way or finite precision arithmetic is 
used. While the first condition is easily ensured in 
practice by using a fixed starting vector for every 
CG inversion, all numerical computations carried 
out using floating point arithmetic are subject to 
rounding errors. 

These rounding errors are normally not con- 
sidered dangerous unless they are exponentially 
amplified. Indeed, without such an amplification, 
the time cost of reducing the error to some preset 
value grows only logarithmically with the num- 
ber of arithmetic operations involved in the com- 
putation. This is a very small correction to the 
growth of the cost of the HMC algorithm as the 
volume and correlation length of the system are 
increased. 

Exponential amplification will occur whenever 
nearby MD trajectories diverge from one another 
exponentially, i.e., when the MD evolution be- 
comes unstable. There are two distinct mecha- 
nisms leading to such an instability ||^. First, this 
is typical for nonlinear equations in the chaotic 
regime. In fact, the existence of a positive lead- 
ing Liapunov exponent for the MD equations of 
pure SU{2) lattice gauge theory was proposed in 
Ref. H] . The second possibility is that the result 
of the discrete integration scheme diverges expo- 
nentially from the true solution. This instability 
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should grow with the number of integration steps 
and is thus expected to have characteristic time 
scale shorter than the one associated with intrin- 
sic chaos. Our numerical results confirm this. 

The integration instability can be analysed in 
the context of free field theory In fact, the be- 
haviour of a single mode with frequency lo already 
reveals all the essential features. One can show 
that the instability accompanied by the exponen- 
tially decaying acceptance rate (P^„ ^ e""^"^) oc- 
cur when ujSt > 2. Here r and St are the tra- 
jectory length and the integration step size re- 
spectively. In Fig. 1 the 'cr — 0' line shows the 
characteristic exponent as a function of St with 
LO fixed to unity. Note the sharp "wall" arising at 
St ^ 2, where the instability sets in. 

Qualitatively similar behaviour is observed for 
the case of many stable modes . The onset of 
instability is determined by the highest frequency 
mode and occurs when uJmaxST = 2. In order to 
keep the acceptance rate constant for free field 
theory as the lattice volume — > cxd, we must de- 
crease St so that VSt'^ stays fixed. Consequently, 
the instabilities go away as we approach the ther- 
modynamic limit. In this sense the leapfrog in- 
stability is a finite volume effect. 

In interacting field theory the notion of inde- 
pendent modes loses its meaning. On the other 
hand, accepting the standard assumption that it 
can be still useful to think in these terms for 
asymptotically free field theories at short dis- 
tances, it is quite plausible to expect similar sce- 
nario there too. The forces acting on the highest 
frequency mode due to the other modes will fluc- 
tuate in some complicated way however, and so 
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Figure 1. Characteristic exponent for the 
fluctuating-frequency harmonic osciUator model. 
The frequency fluctuates around lu = 1. 

we expect that the "wall" at Wi„ax<5T = 2 will get 
smeared out. This is illustrated for the simple 
model of a harmonic oscillator whose frequency 
is randomly chosen from a Gaussian distribution 
with mean to and standard deviation a before 
each MD step. The numerical results shown in 
Fig. 1 confirm that the "wall" in this model does 
indeed spread out. 

Equipped with the above qualitative picture, 
we have studied reversibility numerically for 
SU{3) gauge theory both in the pure gauge and 
dynamical fermion cases Q| (see also related work 
Q). We evolved a typical equilibrium configura- 
tion U using leapfrog equations for some time r, 
then reversed the momenta and evolved it again 
for the same amount of time to get the configura- 
tion U' . Deviations from reversibility were mea- 
sured by 
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but we also recorded the change of energy at the 
end of the trajectory {6H) and at the end of the 
reversed trajectory (ASH). 

In Fig. 2 we collected a typical set of data from 
one pure gauge configuration. The top and bot- 
tom graphs clearly show the integration instabil- 
ity "wall" at St w 0.6, which has spread out as 



expected. At the same time the middle graph in- 
dicates that as we reach the "wall" SH = O(IO^), 
so the integration instabilities are of no practi- 
cal importance for this system. Note however the 
case of the unreasonably long trajectory (r = 40), 
where the reversibility is lost while SH is very 
small implying a good acceptance rate. 
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Figure 2. Results for pure SU{3) gauge theory 
with P — 5.7 on a 4"* lattice as a function of St. 

When ploted as a function of t, all of our 
data show a clear exponential instability in 
1 1 A(5C/ 1 1 . We extracted a characteristic exponent i/ 
(||AJ[/|| e"^) and show the results in Fig. 3. 
Note the same qualitative behaviour we observed 
for the toy model in Fig. 1 except that the in- 
tegration instability "wall" appears at different 
values of St. This probably just reflects the dif- 
ferent highest frequencies of these systems. In 
case of full QCD, the pseudofermions produce a 
force of the order of the inverse lightest fermionic 
mass thus giving the highest relevant frequency 
when simulating close to Kc- This is reflected in 
the bottom graph where the integration instabil- 
ity appears at very small St. 

Notice also that the characteristic exponent 
does not approach zero for small St, which con- 
firms the existence of chaotic continuous time dy- 
namics. Unlike the integration instability, the in- 
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trinsic chaos can not be controlled by adjusting 
6t. Moreover, accepting the standard hypothesis 
that the trajectory length should be scaled pro- 
portionally to the correlation length in order to 
reduce the critical slowing down, non-reversibility 
might cause problems when simulating closer to 
the continuum limit. 
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Figure 3. Characteristic exponent for pure SU{3) 
gauge theory (top), QCD with heavy dynamical 
Wilson quarks (middle), and QCD with light dy- 
namical Wilson quarks (bottom). The data was 
measured on three independent configurations. 

However, our numerical analysis indicates a 
strong /3— dependence of the exponent v, charac- 
terizing the intrinsic chaos. Indeed, Fig. 4 shows 
this for SU{3) pure gauge theory on 4^ and 8* 
lattices. These results can be qualitatively un- 
derstood if we hypothesize that chaos is not only 
a property of this continuous time evolution, but 
is also a property of the underlying continuum 
field theory. This would suggest that i/ scales 
like a physical quantity. At small /? the lattice 
theory is in the strong coupling regime and does 
not obey the asymptotic scaling behaviour. At 
large /3 the system is in a tiny box and is thus 
in the deconfined phase. The finite temperature 
phase transition at Nt = 4 occurs near f3 = 5.7, 
suggesting that the scaling region is in the vicin- 



ity of this value for our lattices. We have fitted 
our 8* data at P — 5.4, 5.5, 5.6, 5.7 to the one 
loop asymptotic scaling form v = ce~^^^^^° , with 
/3q = (11 — |n/)/167r^ and with constant c being 
the only free parameter. 
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Figure 4. Characteristic exponent as a function 
of (3 for pure SU{3) gauge theory. The data was 
measured on three 4* and two 8*^ configurations. 

The resulting fit, shown in Fig. 4 is surprisingly 
good, suggesting that our hypothesis might in- 
deed be correct. This would mean that the char- 
acteristic exponent is constant when measured in 
"physical" units, that is I'S, would be constant 
as ^ oo. If this is the case, then tuning the 
HMC algorithm by varying the trajectory length 
proportionally to the correlation length does not 
lead to any change in the amplification of round- 
ing errors as we simulate closer to the continuum 
limit. 
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